{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Imports\n",
    "\n",
    "# utility modules\n",
    "import glob\n",
    "import os\n",
    "import sys\n",
    "import re\n",
    "\n",
    "# the usual suspects:\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# specialty modules\n",
    "import h5py\n",
    "import pyproj\n",
    "\n",
    "# run matplotlib in 'widget' mode\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import date\n",
    "\n",
    "def diff_2_date(date_1, date_2):\n",
    "    \"\"\"\n",
    "    from w3resource.com\n",
    "    \"\"\"\n",
    "    #date_1\n",
    "    year_1, year_2 = int(date_1[0:4]), int(date_2[0:4])\n",
    "    month_1, month_2 = int(date_1[4:6]), int(date_2[4:6])\n",
    "    day_1, day_2 = int(date_1[6:8]), int(date_2[6:8])\n",
    "    \n",
    "    f_date = date(year_1,month_1, day_1)\n",
    "    l_date = date(year_2,month_2, day_2)\n",
    "    delta = l_date - f_date\n",
    "    return(delta.days)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append(os.path.join(os.getcwd(), '../..'))\n",
    "from readers.read_HDF5_ATL03 import read_HDF5_ATL03\n",
    "from readers.get_ATL03_x_atc import get_ATL03_x_atc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def atl06_to_dict(filename, beam, field_dict=None, index=None, epsg=None):\n",
    "    \"\"\"\n",
    "        Read selected datasets from an ATL06 file\n",
    "\n",
    "        Input arguments:\n",
    "            filename: ATl06 file to read\n",
    "            beam: a string specifying which beam is to be read (ex: gt1l, gt1r, gt2l, etc)\n",
    "            field_dict: A dictinary describing the fields to be read\n",
    "                    keys give the group names to be read, \n",
    "                    entries are lists of datasets within the groups\n",
    "            index: which entries in each field to read\n",
    "            epsg: an EPSG code specifying a projection (see www.epsg.org).  Good choices are:\n",
    "                for Greenland, 3413 (polar stereographic projection, with Greenland along the Y axis)\n",
    "                for Antarctica, 3031 (polar stereographic projection, centered on the Pouth Pole)\n",
    "        Output argument:\n",
    "            D6: dictionary containing ATL06 data.  Each dataset in \n",
    "                dataset_dict has its own entry in D6.  Each dataset \n",
    "                in D6 contains a numpy array containing the \n",
    "                data\n",
    "    \"\"\"\n",
    "    if field_dict is None:\n",
    "        field_dict={None:['latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "    D={}\n",
    "    file_re=re.compile('ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5')\n",
    "    with h5py.File(filename,'r') as h5f:\n",
    "        for key in field_dict:\n",
    "            for ds in field_dict[key]:\n",
    "                if key is not None:\n",
    "                    ds_name=beam+'/land_ice_segments/'+key+'/'+ds\n",
    "                else:\n",
    "                    ds_name=beam+'/land_ice_segments/'+ds\n",
    "                if index is not None:\n",
    "                    D[ds]=np.array(h5f[ds_name][index])\n",
    "                else:\n",
    "                    D[ds]=np.array(h5f[ds_name])\n",
    "                if '_FillValue' in h5f[ds_name].attrs:\n",
    "                    bad_vals=D[ds]==h5f[ds_name].attrs['_FillValue']\n",
    "                    D[ds]=D[ds].astype(float)\n",
    "                    D[ds][bad_vals]=np.NaN\n",
    "    if epsg is not None:\n",
    "        xy=np.array(pyproj.proj.Proj(epsg)(D['longitude'], D['latitude']))\n",
    "        D['x']=xy[0,:].reshape(D['latitude'].shape)\n",
    "        D['y']=xy[1,:].reshape(D['latitude'].shape)\n",
    "    temp=file_re.search(filename)\n",
    "    D['rgt']=int(temp['rgt'])\n",
    "    D['cycle']=int(temp['cycle'])\n",
    "    D['beam']=beam\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190430122344_04920311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20181030210407_04920111_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190730080323_04920411_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190220230230_08320211_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190312235510_11380211_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20181108184743_06280111_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190228224553_09540211_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190623094402_13150311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190506112405_05830311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190611193446_11380311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190620171809_12740311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190823150456_08630411_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190529105941_09340311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190308130329_10700211_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190824143917_08780411_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190812071246_06900411_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20181118033101_07710111_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190822060451_08420411_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190607194306_10770311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /srv/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190129164404_04920211_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "read 613 data files of which 20 gave errors\n"
     ]
    }
   ],
   "source": [
    "# find all the files in the directory:\n",
    "#ATL06_files=glob.glob(os.path.join(data_root, 'PIG_ATL06', '*.h5'))\n",
    "data_root='/srv/shared/surface_velocity/'\n",
    "ATL06_files=glob.glob(os.path.join(data_root, 'FIS_ATL06', '*.h5'))\n",
    "D_dict={}\n",
    "error_count=0\n",
    "for file in ATL06_files:\n",
    "    try:\n",
    "        D_dict[file]=atl06_to_dict(file, '/gt2l', index=slice(0, -1, 25), epsg=3031)\n",
    "    except KeyError as e:\n",
    "        print(f'file {file} encountered error {e}')\n",
    "        error_count += 1\n",
    "print(f\"read {len(D_dict)} data files of which {error_count} gave errors\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "D_2l={}\n",
    "D_2r={}\n",
    "\n",
    "# specify the rgt here:\n",
    "rgt=\"1092\"\n",
    "rgt='0848'\n",
    "# iterate over the repeat cycles\n",
    "for cycle in ['03','04','05','06','07']:\n",
    "    for filename in glob.glob(os.path.join(data_root, 'FIS_ATL06', f'*ATL06_*_{rgt}{cycle}*_003*.h5')):\n",
    "        try:\n",
    "            # read the left-beam data\n",
    "            D_2l[filename]=atl06_to_dict(filename,'/gt1l', index=None, epsg=3031)\n",
    "            # read the right-beam data\n",
    "            D_2r[filename]=atl06_to_dict(filename,'/gt1r', index=None, epsg=3031)\n",
    "            # plot the locations in the previous plot\n",
    "            #map_ax.plot(D_2r[filename]['x'], D_2r[filename]['y'],'k');  \n",
    "            #map_ax.plot(D_2l[filename]['x'], D_2l[filename]['y'],'k');\n",
    "        except Exception as e:\n",
    "            print(f'filename={filename}, exception={e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0e635933303b4e6ca86b1ce3d8329155",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f86ae76ce90>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "file_re=re.compile('ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5')\n",
    "\n",
    "for i, key in enumerate(D_2l.keys()):\n",
    "    if i == 0:\n",
    "        D1 = D_2l[key]\n",
    "        for string in key.split('/'):\n",
    "            if len(string) ==49:\n",
    "                temp=file_re.search(string)\n",
    "        D1['timestamp'] = temp['date']\n",
    "            \n",
    "for i, key in enumerate(D_2l.keys()):\n",
    "    if i == 2:\n",
    "        D2 = D_2l[key]\n",
    "        for string in key.split('/'):\n",
    "            if len(string) ==49:\n",
    "                temp=file_re.search(string)\n",
    "        D2['timestamp'] = temp['date']\n",
    "\n",
    "#for filename, Di in D_2l.items():\n",
    "#Plot only points that have ATL06_quality_summary==0 (good points)\n",
    "ind1=D1['atl06_quality_summary']==0\n",
    "ind2=D2['atl06_quality_summary']==0\n",
    "#define the heights of the segment endpoints.  Leave a row of NaNs so that the endpoints don't get joined\n",
    "h1 = D1['h_li'][ind1]\n",
    "x1 = D1['x_atc'][ind1]\n",
    "h2 = D2['h_li'][ind2]\n",
    "x2 = D2['x_atc'][ind2]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(x1,h1)\n",
    "plt.plot(x2,h2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "96e3ee9bfce541f0b43d990054233b06",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "idx = ~np.isnan(h1) & ~np.isnan(x1) & ~np.isnan(h2) & ~np.isnan(x2)\n",
    "H1 = h1[idx]\n",
    "X1 = x1[idx]\n",
    "Slope1 = np.gradient(H1,X1)\n",
    "H2 = h2[idx]\n",
    "X2 = x2[idx]\n",
    "Slope2 = np.gradient(H2,X2)\n",
    "\n",
    "start, end = 500,1000\n",
    "\n",
    "plt.figure();\n",
    "plt.subplot(211)\n",
    "plt.plot(X1[start:end],H1[start:end], label=\"cycle=3\")\n",
    "plt.plot(X2[start:end],H2[start:end], label=\"cycle=5\")\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('elevation');\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.plot(X1[start:end],Slope1[start:end],ms=1, label=\"cycle=3\")\n",
    "plt.plot(X2[start:end],Slope2[start:end],ms=1, label=\"cycle=5\")\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('slope');\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "90b4a4a3234548c5afd6b214bcd06189",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "velocity of the best fit:  0.04945054945054945 m/d\n",
      "velocity of the best fit:  18.049450549450547 m/a\n"
     ]
    }
   ],
   "source": [
    "idx_offsets = np.arange(-50,50)\n",
    "cos = np.empty_like(idx_offsets).astype(float)\n",
    "for i,idx_offset in enumerate(idx_offsets):\n",
    "    co = np.corrcoef(Slope1[1000:2000],Slope2[1000+idx_offset:2000+idx_offset])\n",
    "    cos[i] = co[1,0]\n",
    "    \n",
    "plt.figure()\n",
    "plt.plot(idx_offsets,cos)\n",
    "\n",
    "idx_best = idx_offsets[np.argmax(cos)]\n",
    "diff_time_days = diff_2_date(D1['timestamp'], D2['timestamp'])\n",
    "\n",
    "print('velocity of the best fit: ' , abs(idx_best)/diff_time_days, 'm/d')\n",
    "print('velocity of the best fit: ' , (abs(idx_best)/diff_time_days)*365, 'm/a')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "34c6381205874b3a9da26802b0967e6b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure();\n",
    "plt.subplot(211)\n",
    "plt.plot(X1[start:end],H1[start:end], label=\"cycle=3\")\n",
    "plt.plot(X2[start:end],H2[start+idx_best:end+idx_best], label=\"cycle=5\")\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('elevation');\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.plot(X1[start:end],Slope1[start:end],ms=1, label=\"cycle=3\")\n",
    "plt.plot(X2[start:end],Slope2[start+idx_best:end+idx_best],ms=1, label=\"cycle=5\")\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('slope');\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
